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Abstract. We here discuss the outcome of an hypothetic experiments of populations 
dynamics, where a set of independent realizations is made available. The importance 
of ensemble average is clarified with reference to the registered time evolution of 
key collective indicators. The problem is here tackled for the logistic case study. 
Theoretical prediction are compared to numerical simulations. 
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1. Introduction 

The problem of explaining the emergence of self-organized, macroscopic, patterns from 
a limited set of rules governing the mutual interaction of a large assembly of microscopic 
actors, is often faced in several domains of physics and biology This challenging task 
defines the realm of complex systems, and calls for novel paradigms to efficiently intersect 
distinct expertise. 

Population dynamics has indeed attracted many scientists pQ and dedicated models 
were put forward to reproduce in silico the change in population over time as displayed 
in real ecosystems (including humans). Two opposite tendencies are in particular to be 
accomodated for. On the one hand, microscopic agents do reproduce themselves with 
a specific rate r, an effect which translates into a growth of the population size P. On 
the other, competition for the available resources (and death) yields a compression of 
the population. In a seminal work by Verhulst [2], these ingredients were formalized in 
the differential equation: 



K is the so called carrying capacity and identifies the maximum allowed population for 
a selected organism, under specific environmental conditions. The above model predicts 
an early exponential growth, which is subsequently antagonized by the quadratic 
contribution, responsible for the asymptotic saturation. The adequacy of the Verhulst's 
model was repeatedly tested versus laboratory experiments: Colonies of bacteria, yeast 
or other simple organic entities were grown, while monitoring the time evolution of the 
population amount. In some cases, an excellent agreement [3j H] with the theory was 
reported, thus supporting the biological validity of Eq. ([T]) . Conversely, the match with 
the theory was definetely less satisfying for e.g. fruit flies, flour beetles and in general for 
other organisms that rely on a more complex life cycle. For those latter, it is necessary 
to invoke a somehow richer modelling scenario which esplicitly includes age structures 
and time delayed effects of overcrowding population [I]. For a more deailed account 
on these issues the interested reader can refer to the review paper [3j and references 
therein. 

Clearly, initial conditions are crucial and need to be accurately determined. An 
error in assessing the initial population, might reflect in the estimates of the parameters 
r and K, which are tuned so to adjust theoretical and experimental data. In general, 
the initial condition relative to one specific experimental realization could be seen 
as randomly extracted from a given distribution. This, somehow natural, viewpoint 
is elaborated in this paper and its implications for the analysis of the experiments 
thoroughly explored. 

In particular we shall focus on the setting where N independent population 
communities are (sequentially or simultaneously) made to evolve. The experiment here 
consists in measuring collective observables, as the average population and associated 
momenta of the ensemble distribution. As anticipated, sensitivity to initial condition do 




(1) 
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play a crucial role and so need to be properly addressed when aiming at establishing a 
link with (averaged) ensemble measurements, or, equivalently, drawing reliable forecast. 
To this end, we will here develop two analytical approaches which enable us to 
reconstruct the sought distribution. The first, to which section [2] is devoted, aims 
at obtaining a complete description of the momenta, as e.g. the mean population 
amount. This is an observable of paramount importance, potentially accessible in real 
experiments. The second, discussed in section HI introduces a master equation which 
rules the evolution of the relevant distribution. It should be remarked that this latter 
approach is a priori more general then the former, as the momenta can in principle be 
calculated on the basis of the recovered distribution. However, computational difficulties 
are often to be faced which make the analysis rather intricate. In this perspective the 
two proposed scenario are to be regarded as highly complementary. 

In the following, for practical purposes, we shall assume each population to evolve 
as prescribed by a Verhulst type of equation. The methods here developed are however 
not limited to this case study but can be straightforwardly generalized to settings were 
other, possibly more complex, dynamical schemes are put forward. 



2. On the momenta evolution 



Imagine to label with Xi the population relative to the i-th realization, belonging to the 
ensemble of N independent replica. As previosuly recalled, we assume each x^ to obey 
a first order differential equation of the logistic type, namely: 

^=x i (l-z i ), (2) 

that can be straightforwardly obtained from (DQ) by setting x = P/K and renaming the 
time t — > rt. The initial condition will be denoted by X®. 

A natural question concerns the expected output of an hypothetic set of experiments 
constrained as above. More concretely, can we describe the distribution of possible 
solutions, once the collection of initial data is entirely specified? 

The m-th momentum associated to the discrete distribution of N repeated 
measurements acquired at time t reads: 

<im>(i) = M)r±__±(^r (3 ) 

To reach our goal, we introduce the time dependent moment generating function, 

oo 

G(£,t):=Y,C <x m >(t). (4) 

m=l 

This is a formal power series whose Taylor coefficients are the momenta of the 
distribution that we are willing to reconstruct, task that can be accomplished using 
the following relation: 
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By exploiting the evolution's law for each Xj, we shall here obtain a partial 
differential equation governing the behavior of G. Knowing G will eventually enables 
us to calculate any sought momentum via multiple differentiation with respect to £ as 
stated in (jSJ). 

Deriving ([3]) and making use of Eq. (T5]) immediately yields: 



— ^ T- m ^ ( + \ — — — — ^-A 

dt >[t) ~ dt ~ N^ X * dt 

i=i i=i 



N 

m 



N 



xT~ l x l (l - Xi) = m (< x m > - < x m+1 >) , (6) 



i=i 



On the other hand, by differentiating (j3j) with respect to time, one obtains : 

§ = E r-^l^ = E ™r (< - m > - < - m+1 >) > ( 7 ) 

m>l m>l 

where used has been made of Eq. ( ffl) . We can now re-order the terms so to express the 
right hand side as a function of G |i] and finally obtain the following non-homogeneous 
linear partial differential equation: 

d t G-(£-l)d 5 G=j. (10) 

Such an equation can be solved for £ close to zero (as in the end of the procedure 
we shall be interested in evaluating the derivatives at £ = 0, see Eq. (jSJ) ) and for all 
positive t. To this end we shall specify the initial datum: 

G(£,0) = ^r < x m > (0) = $(£) , (11) 

m>l 

i.e. the initial momenta or their distribution. 

Before turning to solve (flOl) . we first simplify it by introducing 

G = e 9 namely g = \ogG, (12) 
| Here the following algebraic relations are being used: 

£%G(£, t)=e^ mC- 1 < x m >= ]T mf < x m > , (8) 

rn > 1 m > 1 

and 

^ G^t) = ^ ^ ^ m _j < ^ m >= ^ (m _ <x m > 
^ m>l m>l 

= J2{m- i)e m_1 < x m > 

m>l 

Renaming the summation index, m — 1 — > m, one finally gets (note the sum still begins with m = 1): 

= m£ m < x m+1 > . (9) 

3 ™ \ 1 
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then for any derivative 

9»G = Gd*g , (13) 
where * = £ or * = t, thus ( ITUi) is equivalent to 

% - (£ - 1)% = i , (14) 

with the initial datum 

<?(£,o) = m = log $(£). (is) 

This latter equation can be solved using the method of the characteristics, here 
represented by: 

§ = -«-D, (16) 
which are explicitly integrated to give: 

£(f) = l + (£(0)-l)e-*, (17) 
where £(0) denotes at t = 0. Then the function u(£(t),t) defined by: 

u(m,t) ■= <K£(0)) + J 1 + (e(Q |_ 1)e _ 5 ds , (18) 

is the solution of (TT4"j) . restricted to the characteristics. Observe that u(£(0),0) = 
0(£(O)), so ffTS"]) solves also the initial value problem. 

Finally the solution of ()15p is obtained from u by reversing the relation between 
f(t) and £(0), i.e. £(0) = - l)e* + V. 

= 0((£-l)e* + l)+A(£,t), (19) 

where A(£, t) is the value of the integral in the right hand side of i fTBl . 

This integral can be straightforwardly computed as follows (use the change of 
variable z = e~ s ): 

[ l 1 f e ^ -dz 1 

A = 7o 1 + W) - ^ = Ji ~ T+WTfy ' (20) 

which implies 

a = - r * g - i + S)"-dJ ^ - ,ogz + iog<i + k<o) - i,2) 

= ( + log(l + K(0) - l)e-') - log{(0) . (21) 
According to ( TT91 the solution g is then 

<?(£, t) = <P ((£ - l)e* + 1) + t + log£ - log((f - l)e* + 1) , (22) 
from which G straightforwardly follows: 
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As anticipated, the function G makes it possible to estimate any momentum (j5J). 
As an example, the mean value correspond to setting m — 1, reads: 

<x>{t) = d ( G q = [$' (1 + (£ - l)e*) e* ^ 



+ $(! + (£- l)e*) e' 



(f - l)e* + 1 
(f - l)e* + 1 - £e*n 
(l + (f-l)e*) 2 J 5=0 



e 



7 $(1 - e*) . (24) 
1 — e ( 

In the following section we shall turn to considering a specific application and test 
the adequacy of the proposed scheme. 

3. Uniform distributed initial conditions 

In this section we will focus on a particular case study in the aim of clarifying the 
potential interest of our findings. The inital data (i.e. initial population amount) are 
assumed to span uniformly a bound interval [a, b] . No prior information is hence available 
which favours one specific choice, all accessible initial data being equally probably. To 
fix the ideas we shall here set a = and b = 1/2. The probability distribution ip{x) 
clearly reads 0: 

/n J 2 * * e [0, 1/2] ( . 

W = < _ + , > (25) 

I otherwise 

and cosequently the initial momenta are: 

■1 f l/2 1 , 



< x m > (0) = / mm = / 2C d£ = — -T^ • (26) 
Jo Jo m + i z 

Hence the function $ as defined in ( TTTj) takes the form: 

1 £ m 

$ (£) ~ Yl m+ i2™ ' ^ 

m>l 

A straightforward algebraic manipulation allows us to re-write ( J27l) as follows: 

^ 7/ m I /"J/ ^ i ~ i 

E^i = -/ E- m ^ = -/ r ^^ = -i--io g (i-y),(28) 

m >i m+i ^ m>l ^ 1 Z ^ 

thus 

^(0 = -l-|log (i-f) • (29) 

We can now compute the time dependend moment generating function, G(£,t), 
given by ( 1231) as: 



(e - l)e* + 1 



-1- - log(l + 1 

(f - l)e* + 1 



,(30) 



§ We hereby assume to sample over a large collection of independent replica of the system under 
scrutiny (N is large). Under this hypothesis one can safetly adopt a continuous approximation for 
the distribution of allowed initial data. Conversely, if the number of realizations is small, finite size 
corrections need to be included. 
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and thus recalling (j5J) we get 

< x > (t) = - 



2e* 



< x 2 > (t) 



- 1 

;2t 



-I) 2 

4e 2 * 



log 



e* + 1 



(31) 



log 



2e 



2f 



e* - l) 2 (e* 



1 



(e* - l) 2 (e* - l) 3 ~° V 2 
For large enough times, the distribution of the experiments' outputs is in fact 
concentrated around the asymptotic value 1 with an associated variance (calculated 
from the above momenta) which decreases monotonously with time. In Fig. [1] direct 
numerical simulations are compared to the analytical solution (13Tb ). returning a good 
agreement. A naive approach would suggest interpolating the averaged numerical profile 
with a solution of the logistic model whose initial datum x° acts as a free parameter 
to be adjusted to its best fitted value. As testified by visual inspection of Fig. [1] this 
procedure yields a significant discrepancy, which could be possibly misinterpreted as a 
failure of the underlying logistic evolution law. For this reason, and to avoid drawing 
erroneous conclusions when ensemble averages are computed, attention has to be payed 
on the role of initial conditions. 




time 



Figure 1. Main panel: Time evolution of the first moment < x > ( i ). The (blue) 
solid line stands for direct simulations averaged over N = 100 independent realizations. 
The (green) dashed line represents the analytical solution ([31k ). The (red) dot-dashed 
line is the solution of the logistic Eq. where the initial datum is being adjusted 
to the best fit value x° = 0.2f6. Inset: the solid (resp. dashed) line represents the 
difference between the analytical (resp. fitted) and numerical curves. 



Remark 3.1 (Best parameters estimates). In the preceding discussion the role of initial 
condition was elucidated. In a more general setting one might imagine r, the logistic 
parameter, to be an unknown entry to the model (see Eq. $M)). One could therefore 
imagine to proceed with a fitting strategy which adjusts both x° and r so to match the 
(averaged) data. Alternatively, and provided the distribution of initial conditions is 
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assigned (here assumed uniform), one could involve the explicit solution l[JJ]q) where 
time is scaled back to ist original value: 



-,rt 



< X > (t) 



2e 



ri 



oft 



- 1 



log 



+ 1 



(32) 



( e rt _ !)2 

and let the solely parameter r to run freely so to search for the optimal agreement with 
the data. As an example, we perfomed N = 100 repetead numerical simulations of the 
logistic model with parameter r = 1.5 and intial data uniformly distributed in [0,1/2]. 
Using the straightforward solution of the logistic equation where x° and r are adjusted, 
returns r = 1.2123. The analysis based on (SOfy leads to r = 1.5662, which is definitely 
closer to the true value. 

Remark 3.2 (On the case of a normal distribution). The above discussion is rather 
general and clearly extends beyond the uniform distribution case study. The analysis 
can be in fact adapted to other settings, provided the distribution of initially allowed 
population amount is known. We shall here briefly discuss the rather interesting case 
where a normal distribution is to be considered. Let us assume that x® are random 
normally distributed values with mean 1/4 and standard deviation a 2 , one can compute 
all the intial momenta < x m > (0) as: 

1 -w 



<x m > (0) 



1 f €-1/4 ' 

e 21 



(33) 



0-V27T 

Assuming a , k > 3 to be negligible with respect to a 2 , the function $(£) specifying the 
initial datum in Eq. < f77]j reads: 



m>l m>3 

Collecting together the terms (£/4) m for m > 1 we obtain: 



m(m 



m-2 



2 <-2 



.(34) 



m>l 



(35) 



while the remaining terms read: 

m(m — 1) 



E 

m>3 



4V 



E 

m>2 



m(m — 1) 



4V 



(36) 



It is then easy to verify that their contributution to the required $ funcion results in 



£ , (4a) 2 2(£/4) 2 



4 3 £V 2 



+ 



(37) 



4-^ 2 (l-e/4) 3 4-^ (4-0 3 ' 
To proceed further we again calculate the derivatives of G (defined through the 



function <&), evaluate them at £ 
time, for all m > 1 . 



0, and eventually get the evolution of < x m > in 
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4. Monitoring the time evolution of the probability distribution function of 
expected measurements 

As opposed to the above procedure, one may focus on the distribution function of 
expected outputs, rather then computing its momenta. The starting point of the 
analysis relies on a generalized version of the celebrated Liouville theorem. This latter 
asserts that the phase-space distribution function / is constant along the trajectory 
of the system. For a non Hamiltonian system this condition results in the following 



equation (for convenience derived in the Appendix Appendix A ) for the evolution of the 
probability density function under the action of a generic ordinary differential equation, 
here represented by the vector field X: 

^ + Vf-X + fdivX = 0, (38) 

where divX = Y^dXi/dxi. 

For the case under inspection the 1-dim vector field reads X(x) = x(l — x) and 
hence divX = (1 — 2x). Thus, introducing F = log / Eq. (138j) can be cast in the form: 
OF dF 

— + x(l-x)— + (l-2x) = 0. (39) 

To solve this equation we use once again the methods of characteristics, which are 
now solutions of x = x(l — x), namely: 

/ \ x(Q)e t 

x(t) = -V — r-r^ > 40 

W 1 - x(0) + z(0)e* ' V ; 

The solution of (1391) is hence: 

F(x,t) = F o (x(0)) - f (l-2x(s))ds, (41) 
Jo 

where F = log ip is related to the probability distribution function at t — and must be 
evaluated at x(0), seen as a function of x(t). The integral can be computed as follows: 

ly^^JX 1 - 2 —* ^+ 2 Ml-(0) + *(0)e<) .(42) 
Such an expression has to be introduced into (jUJ) once we explicit x(0) for x(t) = x as: 



xe t 



x(0) = 7 • (43) 

w 1 - x + xe"* v ' 



Hence: 



(xe * \ 
■ -t-21og(l-x + xe-*) , (44) 
1 — x + xe" 1 J v ' 

and finally back to the original /: 



xe 1 \ e 1 



f(x,t)=ip[- 2, (45) 

K ' ^yi-x + xe-tj (\- x + xe -tf v ' 

which stands for the probability density function which describes for all t the expected 
distribution of x's. In Fig.[2]we compare the analytical solutions (I45p with the numerical 
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7 




Figure 2. Time evolution of the probability distribution function. Histograms refers 
to numerical simulation and arc calculated at different time: t = (green online), 
t = i.5 (red, online), t = 2.0 (blue online). The lines represent the corresponding 
analytical solution 

simulation of the logistic model ([2]) under the assumption of N = 1000 initial data 
normally distributed with mean 1/2 and variance 0.005. 

Notice that having calculated the distribution / will enable in turn, at least in 
principle, to to calculate all the associated momenta. 

5. Conclusion 

Forecasting the time evolution of a system which obeys to a specifc governing differential 
equation and is initialized as follows a specific probability distribution, constitutes 
a central problem in several domains of applications. Assume for instance a set of 
independent measurements to return an ensemble average which is to be characterized 
according to a prescribed model. Biased conclusion might result from straightforward 
fitting strategies which do not correctly weight the allowed distribution of initial 
condition. 

In this paper we address this problem by providing an exact formula for the time 
evolution of momenta and probability distribution function of expected measurements, 
which is to be invoked for a repeaded set of indipendent experiments. Though general, 
the method is here discussed with reference to a simple, demonstrative problem of 
population dynamics. 
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Appendix A. The generalized Liouville theorem 

Let X(x) be a vector field to which we associate the ordinary differential equation: 

x = X(x) VxGfi, (A.l) 

where Q is the phase space. Suppose to define a probability density function of the 
initial data on Q. Namely we have a function ip defined in the phase space Q, such that 
for all Bcfl, f B ip(x)dx denotes the probability that a randomly drawn initial datum 
will belong to B and f n ip(x)dx = 1. 

We are interested in determining for any t > 0, the probability that a solution 
of flA.lj) will fall in a open set B' C Q. Let us call f(x, t) such probability, by continuity 



we must have f(x, 0) = ip(x) and J n f(x, t)dx = 1 for all t > 0. 

For any B C Q, P(B) = f B f(x,t)dx denotes the probability to find a point in B 
at time t. We can then assume that this probability does not change if the set B' is 
transported by the flow of (1A.1I) . P(B) = P(A) where A = $ S (-B), being $ s the flow at 
time s of the vector field. Namely 



f(y,t + s)dy = / f(x,t)dx, (A.2) 

A=<l> a (B) JB 

the change of coordinates y = $ s (x) allows to rewrite the previous relation as follows: 



f(y,t+s)dy= / f($ s (x),t+s) det D$ s (x)dx = / f(x, t) dx , (A.3) 

A=<S> S {B) JB JB 

being DQ s (x) the Jacobian of the change of variables. 

The relation (1A.3I) should be valid for any set B, thus: 

f(x, t) = /($ s (x), t + s) det D$ s (x) , (A.4) 

for all x G Q and for all t, s. Deriving with respect to s and evaluating the derivative at 
s = we get the required relation (recall D$°(x) = identity): 

^(x, t) + VJ(x, t) ■ X(x) + f(x, t)divX(x) = . (A.5) 
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